% -- A Ramsey Theory of Financial Distortions --
% Plotting Figures 1 and 2
clear; close all; clc

%% ----------------------  Parameters   -----------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% NOTE: many parameters below allow generalization of the two-period
%%%% model. The values (such as full depreciation par.DELTA = 1) correspond 
%%%% to the simple model in the paper.
par.BETA    = 0.96;
par.phi     = 0.5;
par.ALPHA   = 0.33;
par.DELTA   = 1; % depreciation rate
par.MU      = 1; % disutility of labor supply
par.NU      = 1; % inverse frisch
par.A       = 1; % aggregate productivity
par.G_1     = 0; % government spending in period 1
par.G_2     = 0; % government spending in period 2


%% ----------------------  Preparation  -----------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PSI_start = 0;
PSI_end   = 1;
K_c_0   = (par.A * (1 - par.ALPHA) / par.MU) ^(1 / par.NU) * ((1 / par.BETA - (1 - par.DELTA)) / par.A / par.ALPHA)^((par.ALPHA + par.NU) / (par.ALPHA - 1) / par.NU);
K_c_inf = (par.A * (1 - par.ALPHA) / (1 + par.NU) / par.MU)^(1 / par.NU)...
    * ((1 / par.BETA / par.phi - (1 - par.DELTA)) / par.A / par.ALPHA)^((par.ALPHA + par.NU) / (par.ALPHA - 1)/ par.NU);
K_u_inf = (par.A * (1 - par.ALPHA) / (1 + par.NU) / par.MU)^(1 / par.NU)...
    * ((1 / par.BETA - (1 - par.DELTA)) / par.A / par.ALPHA)^((par.ALPHA + par.NU) / (par.ALPHA - 1)/ par.NU);

%%%% bounds for B_e_0

display('The highest possible liquidity level so that the economy is always financially constrained')
B_e_0_bound_1 = (1 - par.phi) * K_c_inf 

display('The highest possible liquidity level so that the economy is never unconstrained')
B_e_0_bound_2 = (1 - par.phi) * K_u_inf 

display('The lowest possible liquidity level so that the economy is never constrained')
B_e_0_bound_3 = (1 - par.phi) * K_c_0


%% ----------------------  Figure 1  ---------------------------------------
display(' ---------   Figure 1: The constrained and unconstrained solution   ---------------')

B_e_0_example = 0.026; %%%% example liquidity level that features constrained
K_star_example = B_e_0_example / (1 - par.phi);

display('The threshold that divides the constrained and the kink regions')
PSI_c = fsolve(@(PSI)cut_off_PSI(PSI, B_e_0_example, par, 'constrained'), 0.1)
display('The threshold that divides the constrained and the unconstrained regions')
PSI_u = fsolve(@(PSI)cut_off_PSI(PSI, B_e_0_example, par, 'unconstrained'), 0.2)

figure(1)
fplot(@(PSI) ((1 + PSI) ./ (1 + (1 + par.NU) * PSI) * par.A * (1 - par.ALPHA) / par.MU).^(1 / par.NU) ...
    .* ( ((par.phi + PSI) / par.phi ./ (1 + PSI) / par.BETA - (1 - par.DELTA)) / par.A / par.ALPHA) .^((par.ALPHA + par.NU) / (par.ALPHA - 1) / par.NU), [PSI_start, PSI_c], '-', 'LineWidth', 5, 'color', 'r'), hold on;

fplot(@(PSI) ((1 + PSI) ./ (1 + (1 + par.NU) * PSI) * par.A * (1 - par.ALPHA) / par.MU).^(1 / par.NU) ...
    .* ( ((par.phi + PSI) / par.phi ./ (1 + PSI) / par.BETA - (1 - par.DELTA)) / par.A / par.ALPHA) .^((par.ALPHA + par.NU) / (par.ALPHA - 1) / par.NU), [PSI_start, PSI_end], '--', 'LineWidth', 2, 'color', 'b')

xlabel('\Psi_1','FontSize', 16)
ylim([0.02, 0.07])
txt = {'K^c(\Psi_1)'};
text((PSI_start + PSI_end) / 2, 0.03,txt, 'FontSize',16)

fplot(@(PSI) ((1 + PSI) ./ (1 + (1 + par.NU) * PSI) * par.A * (1 - par.ALPHA) / par.MU).^(1 / par.NU) ...
   * ( (1 / par.BETA - (1 - par.DELTA)) / par.A / par.ALPHA) .^((par.ALPHA + par.NU) / (par.ALPHA - 1) / par.NU), [PSI_u, PSI_end], 'LineWidth', 5, 'color', 'r'); hold on
fplot(@(PSI) ((1 + PSI) ./ (1 + (1 + par.NU) * PSI) * par.A * (1 - par.ALPHA) / par.MU).^(1 / par.NU) ...
    * ( (1 / par.BETA - (1 - par.DELTA)) / par.A / par.ALPHA) .^((par.ALPHA + par.NU) / (par.ALPHA - 1) / par.NU), [PSI_start, PSI_end], '-.', 'LineWidth', 2, 'color', 'b'); hold on
txt = {'K^u(\Psi_1)'};
text((PSI_start + PSI_end) / 2, 0.045,txt,'FontSize',16)

yline(K_star_example, ':', 'K^*', 'Fontsize', 16, 'LineWidth',1.5);
fplot(@(PSI) PSI - PSI + K_star_example, [PSI_c, PSI_u], 'color', 'r', 'LineWidth', 5)

text(0.10,0.022,'\Psi^c','FontSize',14)
text(0.47,0.022,'\Psi^u','FontSize',14)

%% ----------------------  Figure 2  ---------------------------------------

display(' ---------   Figure 2: Plot the Ramsey solution for different B_e_0  ---------------')

%%%% Scenario 1: always financially constrained
grid_number = 40;

B_e_0_1 = B_e_0_bound_1;
PSI_step_1 = 0.05;
[K_1_1, tau_k_2_1, L_1_1, L_2_1, q_1_1, B_0_1, tau_l_2_1] = eqm_constrained(B_e_0_1, par, grid_number, PSI_step_1);
PSI_label_1 = 0:PSI_step_1:PSI_step_1 * (grid_number - 1);

%%%% Scenario 2: financially constrained for small PSI, then kink, and then at the unconstrained region
display('Using a particular value for liquidity held by the entreprenuers')
B_e_0_2 = (B_e_0_bound_3 + B_e_0_bound_2) / 2 
PSI_star_2_constrained      = fsolve(@(PSI)cut_off_PSI(PSI, B_e_0_2, par, 'constrained'),1);
PSI_star_2_unconstrained    = fsolve(@(PSI)cut_off_PSI(PSI, B_e_0_2, par, 'unconstrained'),1);

grid_number_constrained = 101;
grid_number_kink = 20;
grid_number_unconstrained = 150;

PSI_step_2 = PSI_star_2_constrained / grid_number_constrained;
[K_1_2_con, tau_k_2_2_con, L_1_2_con, L_2_2_con, q_1_2_con, B_0_2_con, tau_l_2_2_con] = eqm_constrained(B_e_0_2, par, grid_number_constrained, PSI_step_2);
PSI_label_2_constrained = 0:PSI_step_2:PSI_step_2 * (grid_number_constrained - 1);    

PSI_step_2 = (PSI_star_2_unconstrained - PSI_star_2_constrained) / grid_number_kink;
[K_1_2_kink, tau_k_2_2_kink, L_1_2_kink, L_2_2_kink, q_1_2_kink, B_0_2_kink, tau_l_2_2_kink] = eqm_kink(B_e_0_2, par, grid_number_kink, PSI_step_2, PSI_star_2_constrained);
PSI_label_2_kink = PSI_star_2_constrained:PSI_step_2:(PSI_star_2_constrained + PSI_step_2 * (grid_number_kink - 1));

PSI_step_2 = 0.01;
[K_1_2_uncon, tau_k_2_2_uncon, L_1_2_uncon, L_2_2_uncon, q_1_2_uncon, B_0_2_uncon, tau_l_2_2_uncon] = eqm_unconstrained(B_e_0_2, par, grid_number_unconstrained, PSI_step_2, PSI_star_2_unconstrained);                                                                   
PSI_label_2_unconstrained = PSI_star_2_unconstrained:PSI_step_2:(PSI_star_2_unconstrained + PSI_step_2 * (grid_number_unconstrained - 1));

K_1_2 = [K_1_2_con;K_1_2_kink;K_1_2_uncon]; 
tau_k_2_2 = [tau_k_2_2_con; tau_k_2_2_kink; tau_k_2_2_uncon]; 
tau_l_2_2 = [tau_l_2_2_con; tau_l_2_2_kink; tau_l_2_2_uncon]; 
L_1_2 = [L_1_2_con; L_1_2_kink; L_1_2_uncon]; 
L_2_2 = [L_2_2_con; L_2_2_kink; L_2_2_uncon];
q_1_2 = [q_1_2_con; q_1_2_kink; q_1_2_uncon];
B_0_2 = [B_0_2_con; B_0_2_kink; B_0_2_uncon];
PSI_label_2 = [PSI_label_2_constrained, PSI_label_2_kink, PSI_label_2_unconstrained];

%%%%
figure(2)
subplot(2,3,1)
h = plot(PSI_label_1, K_1_1, PSI_label_2, K_1_2); grid on;
set(h(1),'linewidth',2,'linestyle','-')
set(h(2),'linewidth',2,'linestyle','--')
ylim([0 max(K_1_1)])
set(get(gca,'ylabel'),'Rotation',0)
set(gca,'FontSize', 22)
title('$K_1$','Fontsize',  18, 'interpreter','latex')

subplot(2,3,2)
h = plot(PSI_label_1, L_2_1, PSI_label_2, L_2_2); grid on;
set(h(1),'linewidth',2,'linestyle','-')
set(h(2),'linewidth',2,'linestyle','--')
set(gca,'FontSize', 22)
ylim([0 max(L_2_1)])
title('$L_2$','Fontsize',  18, 'interpreter','latex')

subplot(2,3,3)
h = plot(PSI_label_1, q_1_1, PSI_label_2, q_1_2); grid on;
set(h(1),'linewidth',2,'linestyle','-')
set(h(2),'linewidth',2,'linestyle','--')
set(gca,'FontSize', 22)
title('$q_1$', 'Fontsize',  18, 'interpreter','latex')

subplot(2,3,4)
h = plot(PSI_label_1, tau_k_2_1 * 100, PSI_label_2, tau_k_2_2 * 100); grid on;
set(h(1),'linewidth',2,'linestyle','-')
set(h(2),'linewidth',2,'linestyle','--')
xlabel('$\Psi_1$','Fontsize',  22, 'interpreter','latex')
ylabel('%')
ylim([-90 20])
set(gca,'FontSize', 22)
title('$\tau^{k}_2$', 'Fontsize', 18, 'interpreter','latex')

subplot(2,3,5)
h = plot(PSI_label_1, tau_l_2_1 * 100, PSI_label_2, tau_l_2_2 * 100); grid on;
set(h(1),'linewidth',2,'linestyle','-')
set(h(2),'linewidth',2,'linestyle','--')
xlabel('$\Psi_1$','Fontsize', 22, 'interpreter','latex')
ylabel('%')
set(gca,'FontSize', 22)
title('$\tau^{\ell}_2$', 'Fontsize', 18, 'interpreter','latex')

subplot(2,3,6)
h = plot( PSI_label_1, B_0_1, PSI_label_2, B_0_2); grid on;
set(h(1),'linewidth',2,'linestyle','-')
set(h(2),'linewidth',2,'linestyle','--')
xlabel('$\Psi_1$','Fontsize', 22, 'interpreter','latex')
set(gca,'FontSize', 22)
case_1 = "B^e_0 = "+ compose("%.4f", B_e_0_1);
case_2 = "B^e_0 = "+ compose("%.4f", B_e_0_2);
legend(case_1, case_2, 'location','southeast','Orientation','horizontal','FontSize',14, 'Box','off')
title('$B_0$', 'Fontsize',  18, 'interpreter','latex')
